##
rm(list = ls())
library(ctgt)
library(faux)


#full ct
# get all supsets of R
getsupsets = function(lowset, highset){
  setd = setdiff(highset,lowset)
  setdl = length(setd)
  if(setdl==0){
    return()
  }else{
    allsub = unlist(lapply(setdl:0, combn,  x = setd, simplify = FALSE), recursive = FALSE)
    allsup = lapply(allsub, function(i) c(lowset, i))
    return(allsup)
  }
}

fullct <- function(y, X, xs, hyps){
  supsets = getsupsets(hyps, xs)
  for(i in 1:length(supsets)){
    pv = gt2(y,X,supsets[[i]])[1]
    if(pv>0.05){
      return("not reject")
      break
    }
  }
  return("reject")
}


## size
n = 100

# binary response
y = rbinom(n,1,0.6)


## big design matrix with max 200 features
m = 38
varnames = paste("x", sep = "", 1:m)
nrep = 10
Ctime = matrix(NA,nrep,m/2)
for(s in 1:nrep){
  set.seed(1234 + s)
  X = as.matrix(rnorm_multi(n = n, vars = m, mu = 0, sd = 1, r = 0.5, varnames = varnames) )
  ## first half features are true
  X[which(y==1),1:(m/2)] = X[which(y==1),1:(m/2)] + 0.5
  
  ## features
  vset = seq.int(2,m,2)
  v = length(vset)
  CT = rep(0,v)
  for( i in 1:v){
    hyps = varnames[1:i]
    xs = varnames[c(1:i, (m/2+1):(m/2 +i) ) ]
    ctt = c()
    for(r in 1:10){
      ct_time = system.time( { fct = fullct(y = y, X = X[,xs], xs = xs, hyps = hyps)} )
      ctt[r] = ct_time[1]
      if(ct_time[1] > 300) break
    }
    CT[i] = median(ctt)
    if(CT[i] > 3600) break
  }
  
  Ctime[s,] =  CT
}

Ctime
save(Ctime, file = "CT.RData")